##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--Gene"), type = "character") ,
    make_option(c("--stad_mut_rate_dgc"), type = "character") ,
    make_option(c("--stad_mut_rate_igc"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    Gene <- "TP53"
    work_dir <- ""
    stad_mut_rate_dgc <- paste(work_dir,"/TCGA_STAD_DGC.mut_rate.tsv",sep="")
    stad_mut_rate_igc <- paste(work_dir,"/TCGA_STAD_IGC.mut_rate.tsv",sep="")
	images_path <- paste(work_dir,"/mutRatePlot",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

Gene <- opt$Gene
stad_mut_rate_dgc <- opt$stad_mut_rate_dgc
stad_mut_rate_igc <- opt$stad_mut_rate_igc
images_path <- opt$images_path

###########################################################################################
stad_mut_rate_dgc <- data.frame(fread(stad_mut_rate_dgc))
stad_mut_rate_dgc$SampleNum <- 60
stad_mut_rate_dgc$Class <- "DGC"
stad_mut_rate_igc <- data.frame(fread(stad_mut_rate_igc))
stad_mut_rate_igc$SampleNum <- 150
stad_mut_rate_igc$Class <- "IGC"
stad_mut_rate_dgc <- subset(stad_mut_rate_dgc , gene == Gene)
stad_mut_rate_igc <- subset(stad_mut_rate_igc , gene == Gene)
mut_rate_gene_file <- rbind(stad_mut_rate_dgc,stad_mut_rate_igc)

col <- c(
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255)
  )

names(col) <- c("IGC" , "DGC")

###########################################################################################

dat_plot <- mut_rate_gene_file
dat_plot$Class <- factor( dat_plot$Class , levels = c("IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$rate , 2) * 100 , "%") 

###########################################################################################
## 计算P值
dat_plot$p.value = ""
dat_plot$p_text = ""

for(geneN in unique(dat_plot$gene)){

    print(geneN)

    tmp_1 <- subset( dat_plot , gene == geneN & Class %in% c("IGC") )
    tmp_2 <- subset( dat_plot , gene == geneN & Class %in% c("DGC") )

    if(nrow(tmp_1)==0){
        tmp_1 <- tmp_2
        tmp_1$MutNum <- 0
        tmp_1$MutRate <- 0
        tmp_1$value_text <- ""
    }

    if(nrow(tmp_2)==0){
        tmp_2 <- tmp_1
        tmp_2$MutNum <- 0
        tmp_2$MutRate <- 0
        tmp_2$value_text <- ""
    }

    tmp <- rbind( tmp_1 , tmp_2 )

    tmp_fisher <- matrix(c(tmp$freq , tmp$SampleNum - tmp$freq) , ncol = 2)
    p <- fisher.test(tmp_fisher)$p.value

    dat_plot[dat_plot$gene == geneN & dat_plot$Class %in% c("IGC" , "DGC"),"p.value"] <- p
    p_text <- paste0( "DGC vs IGC \np = " , format(as.numeric(p) , scientific = TRUE , digits = 3) )
    dat_plot[dat_plot$gene == geneN & dat_plot$Class %in% c("IGC" , "DGC"),"p_text"] <- p_text
}

dat_plot$gene <- factor( dat_plot$gene , 
    levels = unique(dat_plot[order(dat_plot$freq , decreasing=T),"gene"]) , order = T)

###########################################################################################

dat_plot <- subset( dat_plot , gene == Gene )

dat_plot2 <- dat_plot
dat_plot2$Type <- ""
dat_plot2$Type[dat_plot2$Class %in% c("IGC" , "DGC")] <- "IGC vs DGC"
dat_plot2$p_text <- paste0( "p = " , format(as.numeric(dat_plot2$p.value) , scientific = TRUE , digits = 3))
dat_plot2$Type <- factor( dat_plot2$Type , levels = "IGC vs DGC" )

###########################################################################################
## IGC vs DGC
plot <- ggplot( data = dat_plot2[dat_plot2$Type=="IGC vs DGC",] , aes( x = Class , y = rate , fill = Class ))+
geom_bar(position = "stack", stat = "identity") + 
theme_bw()+
labs(x="",y="Mutation Rate")+
facet_wrap(~Type,scales="free")+
theme(panel.grid = element_blank())+
scale_fill_manual(values=col) +
ylim(0,1.05)+
geom_text(aes(label=p_text , y = 1 ,x = 1.5),size=3)+
geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , fontface='bold' , color="black")+
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='right',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 7,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            strip.text.x = element_text(size = 7 , face = 'bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 8,color="black",face='bold') ,
            )

out_name <- paste0( images_path , "/TCGA_MutRate_" , Gene , ".IGC_DGC.pdf" )
ggsave( out_name , plot , width = 3.5 , height = 5 )


